Load and install R packages:
packages <- c('lmerTest', 'lme4', 'ggplot2', 'tidyverse', 'readxl', 'purrr', 'performance', 'emmeans', 'MASS')
# Check if each package is installed; if not, install it
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
library(pkg, character.only = TRUE)
}
Grab neophobia_data.csv from processed_data
directory:
script_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)
# Set the working directory to the script's directory
setwd(script_dir)
data <- read.csv("processed_data/neophobia_data.csv", row.names = 1)
Summary statistics of the data:
# Quick look at structure
str(data)
'data.frame': 536 obs. of 15 variables:
$ Bird_ID : chr "BB_BB" "BB_BB" "BB_BB" "BB_BB" ...
$ Trial : int 1 2 3 4 5 6 7 8 1 2 ...
$ Enclosure : chr "B8" "B8" "B8" "B8" ...
$ Object : chr "novel" "novel" "control" "control" ...
$ Context : chr "group" "individual" "group" "individual" ...
$ Latency_to_enter: num 2.2 6.23 2.27 1.57 1.07 ...
$ Latency_to_Eat : num 2.9 13.83 2.2 3.73 2.73 ...
$ Zoi_duration : num 183.7 91.9 85.4 203.3 371.5 ...
$ Object_Type : int 2 4 3 3 1 5 3 3 5 1 ...
$ GroupID : chr "7A" NA "7A" NA ...
$ NestID : chr "2024VM08" "2024VM08" "2024VM08" "2024VM08" ...
$ group_dummy : int 1 0 1 0 1 0 1 0 1 0 ...
$ ind_dummy : int 0 1 0 1 0 1 0 1 0 1 ...
$ Object_contrast : num 0.5 0.5 -0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5 ...
$ Context_contrast: num 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 ...
summary(data)
Bird_ID Trial Enclosure Object Context Latency_to_enter Latency_to_Eat
Length:536 Min. :1.00 Length:536 Length:536 Length:536 Min. : 0.800 Min. : 1.792
Class :character 1st Qu.:2.75 Class :character Class :character Class :character 1st Qu.: 1.467 1st Qu.: 2.800
Mode :character Median :4.50 Mode :character Mode :character Mode :character Median : 1.880 Median : 3.633
Mean :4.50 Mean : 8.380 Mean : 42.457
3rd Qu.:6.25 3rd Qu.: 2.600 3rd Qu.: 5.684
Max. :8.00 Max. :600.000 Max. :600.000
Zoi_duration Object_Type GroupID NestID group_dummy ind_dummy Object_contrast Context_contrast
Min. : 0.00 Min. :1.000 Length:536 Length:536 Min. :0.0 Min. :0.0 Min. :-0.5 Min. :-0.5
1st Qu.: 38.93 1st Qu.:2.000 Class :character Class :character 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:-0.5 1st Qu.:-0.5
Median : 76.35 Median :3.000 Mode :character Mode :character Median :0.5 Median :0.5 Median : 0.0 Median : 0.0
Mean :128.20 Mean :2.787 Mean :0.5 Mean :0.5 Mean : 0.0 Mean : 0.0
3rd Qu.:162.03 3rd Qu.:4.000 3rd Qu.:1.0 3rd Qu.:1.0 3rd Qu.: 0.5 3rd Qu.: 0.5
Max. :598.13 Max. :5.000 Max. :1.0 Max. :1.0 Max. : 0.5 Max. : 0.5
The GroupID values are currently coded as
NA during individual trials. Assign the most frequent group
for each bird:
most_frequent_group <- function(group_ids) {
group_ids <- group_ids[!is.na(group_ids)]
if (length(group_ids) == 0) return(NA)
names(sort(table(group_ids), decreasing = TRUE))[1]
}
most_frequent_group_per_bird <- data %>%
group_by(Bird_ID) %>%
summarize(most_frequent_group = most_frequent_group(GroupID)) %>%
ungroup()
data <- data %>%
left_join(most_frequent_group_per_bird, by = "Bird_ID") %>%
mutate(GroupID = ifelse(is.na(GroupID), most_frequent_group, GroupID)) %>%
dplyr::select(-most_frequent_group)
Adjust the Trial variable by subtracting 1:
data$Trial <- data$Trial - 1
The full model as described in the RR, includes
Object_contrast, Context_contrast, and
Trial, and a complex random effect structure:
enter_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
boundary (singular) fit: see help('isSingular')
summary(enter_model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data
REML criterion at convergence: 5712.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.4529 -0.1408 -0.0422 0.0686 10.2384
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 1.198e+03 34.6124
group_dummy 2.852e-02 0.1689 1.00
NestID (Intercept) 0.000e+00 0.0000
GroupID group_dummy 0.000e+00 0.0000
Residual 2.268e+03 47.6258
Number of obs: 536, groups: Bird_ID, 67; NestID, 43; GroupID, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17.4867 4.3340 279.0124 4.035 7.06e-05 ***
Object_contrast 5.1073 4.1143 464.9955 1.241 0.21510
Context_contrast -12.0511 5.8852 79.6165 -2.048 0.04389 *
Trial -2.6020 0.9052 471.1684 -2.874 0.00423 **
Object_contrast:Context_contrast -9.7670 8.2292 464.9964 -1.187 0.23588
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_ Trial
Objct_cntrs 0.005
Cntxt_cntrs -0.356 0.000
Trial -0.731 -0.007 0.007
Objct_cn:C_ -0.010 0.000 0.000 0.013
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
check_model(enter_model)
The random effect structure seems to be too complex for the amount of
data. However, the variance for the (Intercept) under
NestID is 0 and the variance for
group_dummy under GroupID is 0.
This suggests that both nest as differences between groups contribute
minimally to the variance in the outcome model. Let’s drop both
effects.
enter_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
boundary (singular) fit: see help('isSingular')
summary(enter_model2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data
REML criterion at convergence: 5712.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.4529 -0.1408 -0.0422 0.0686 10.2384
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 1198.0356 34.6127
group_dummy 0.0286 0.1691 1.00
Residual 2268.2114 47.6257
Number of obs: 536, groups: Bird_ID, 67
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17.4867 4.3340 279.0054 4.035 7.06e-05 ***
Object_contrast 5.1073 4.1143 464.9960 1.241 0.21510
Context_contrast -12.0511 5.8852 79.6162 -2.048 0.04389 *
Trial -2.6020 0.9052 471.1689 -2.874 0.00423 **
Object_contrast:Context_contrast -9.7670 8.2292 464.9969 -1.187 0.23588
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_ Trial
Objct_cntrs 0.005
Cntxt_cntrs -0.356 0.000
Trial -0.731 -0.007 0.007
Objct_cn:C_ -0.010 0.000 0.000 0.013
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
check_model(enter_model2)
Seems like the interaction between object and context is non-significant, let’s drop it. Given the non-normal distribution, let’s logtransform the data as well:
enter_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
boundary (singular) fit: see help('isSingular')
summary(enter_model3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data
REML criterion at convergence: 1123.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.9335 -0.4596 -0.0811 0.3538 6.9301
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 0.307211 0.55427
group_dummy 0.008857 0.09411 1.00
Residual 0.386408 0.62162
Number of obs: 536, groups: Bird_ID, 67
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.29304 0.06319 188.76075 20.464 <2e-16 ***
Object_contrast 0.10177 0.05370 465.99922 1.895 0.0587 .
Context_contrast -0.19502 0.07774 76.68931 -2.508 0.0142 *
Trial -0.13771 0.01179 468.78462 -11.679 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_
Objct_cntrs 0.004
Cntxt_cntrs -0.458 0.000
Trial -0.653 -0.007 0.007
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
check_model(enter_model3)
There are still issues with the data distribution, let’s remove birds that did not interact and boxcox_transform:
data_enter <- data %>% filter(Latency_to_enter != 600)
boxcox_transform <- boxcox(lm(Latency_to_enter ~ Object_contrast + Context_contrast + Trial, data = data_enter))
best_lambda <- boxcox_transform$x[which.max(boxcox_transform$y)]
data_enter$Latency_to_enter_trans <- (data_enter$Latency_to_enter^best_lambda - 1) / best_lambda
enter_model4_boxcox <- lmer(Latency_to_enter_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
summary(enter_model4_boxcox)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data_enter
REML criterion at convergence: -115.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.00860 -0.62234 0.05268 0.62220 2.42212
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 0.023257 0.15250
group_dummy 0.007021 0.08379 0.84
Residual 0.037006 0.19237
Number of obs: 533, groups: Bird_ID, 67
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.686462 0.020720 165.959664 33.130 <2e-16 ***
Object_contrast 0.007148 0.016674 395.856100 0.429 0.668
Context_contrast -0.028455 0.020248 64.608541 -1.405 0.165
Trial -0.051896 0.003669 424.210940 -14.146 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_
Objct_cntrs 0.007
Cntxt_cntrs -0.301 -0.002
Trial -0.623 -0.009 0.016
check_model(enter_model4_boxcox)
Comparing the different models, the boxcox transfromed comes out as best:
enter_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
boundary (singular) fit: see help('isSingular')
enter_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
boundary (singular) fit: see help('isSingular')
enter_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
boundary (singular) fit: see help('isSingular')
enter_model4_boxcox <- lmer(Latency_to_enter_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
anova(enter_model, enter_model2, enter_model3, enter_model4_boxcox)
refitting model(s) with ML (instead of REML)
Data: data_enter
Models:
enter_model3: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
enter_model4_boxcox: Latency_to_enter_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
enter_model2: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
enter_model: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
enter_model3 8 959.8 994.0 -471.91 943.8
enter_model4_boxcox 8 -127.5 -93.2 71.73 -143.5 1087.3 0
enter_model2 9 5148.8 5187.3 -2565.41 5130.8 0.0 1 1
enter_model 11 5152.8 5199.9 -2565.41 5130.8 0.0 2 1
A similar model is built for Latency_to_Eat, again
incorporating interaction terms and random effects:
eat_model <- lmer(log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(- 1 + group_dummy | GroupID) +
(- 1 + ind_dummy + group_dummy | Bird_ID),
data = data)
boundary (singular) fit: see help('isSingular')
summary(eat_model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data
REML criterion at convergence: 1636.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.1780 -0.5149 -0.1108 0.1927 4.6862
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 0.78398 0.8854
group_dummy 0.08989 0.2998 1.00
NestID (Intercept) 0.00000 0.0000
GroupID group_dummy 0.00000 0.0000
Residual 1.00669 1.0033
Number of obs: 536, groups: Bird_ID, 67; NestID, 43; GroupID, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.89202 0.10740 165.73052 17.616 < 2e-16 ***
Object_contrast 0.74392 0.08668 465.00491 8.583 < 2e-16 ***
Context_contrast -1.00808 0.11239 90.13817 -8.969 3.94e-14 ***
Trial -0.04274 0.01899 466.49134 -2.251 0.0248 *
Object_contrast:Context_contrast -0.85108 0.17337 465.00516 -4.909 1.27e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_ Trial
Objct_cntrs 0.004
Cntxt_cntrs -0.434 0.000
Trial -0.619 -0.007 0.008
Objct_cn:C_ -0.008 0.000 0.000 0.013
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
check_model(eat_model)
Similarly, the variance for the (Intercept) under
NestID is 0 and the variance for
group_dummy under GroupID is 0.
Let’s drop those:
eat_model2 <- lmer(log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
(- 1 + ind_dummy + group_dummy | Bird_ID),
data = data)
summary(eat_model2)
check_model(eat_model2)
There are still some issues with the normality. Let’s remove non participating birds, and try the same boxcox tranformation:
boxcox model seems to fit data best
eat_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
boundary (singular) fit: see help('isSingular')
eat_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
boundary (singular) fit: see help('isSingular')
eat_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
boundary (singular) fit: see help('isSingular')
eat_model4_boxcox <- lmer(Latency_to_eat_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
anova(eat_model, eat_model2, eat_model3, eat_model4_boxcox)
refitting model(s) with ML (instead of REML)
Data: data_eat
Models:
eat_model3: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
eat_model4_boxcox: Latency_to_eat_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
eat_model2: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
eat_model: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
eat_model3 8 690.1 723.9 -337.03 674.1
eat_model4_boxcox 8 -794.8 -761.0 405.39 -810.8 1484.8 0
eat_model2 9 4137.0 4175.0 -2059.50 4119.0 0.0 1 1
eat_model 11 4141.0 4187.5 -2059.50 4119.0 0.0 2 1
The full model as described in the RR, includes
Object_contrast, Context_contrast, and
Trial, and a complex random effect structure:
summary(enter_model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data_enter
REML criterion at convergence: 4110
Scaled residuals:
Min 1Q Median 3Q Max
-5.7468 -0.0958 -0.0198 0.0611 19.4975
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 1.196e+02 10.93520
group_dummy 3.693e-04 0.01922 1.00
NestID (Intercept) 0.000e+00 0.00000
GroupID group_dummy 0.000e+00 0.00000
Residual 1.722e+02 13.12072
Number of obs: 506, groups: Bird_ID, 67; NestID, 43; GroupID, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.2873 1.2911 247.2473 4.095 5.72e-05 ***
Object_contrast -0.8899 1.1804 435.9808 -0.754 0.4513
Context_contrast -1.7996 1.7837 70.9317 -1.009 0.3164
Trial -0.5942 0.2595 438.0996 -2.289 0.0225 *
Object_contrast:Context_contrast 1.9819 2.3611 435.9820 0.839 0.4017
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_ Trial
Objct_cntrs 0.057
Cntxt_cntrs -0.433 -0.040
Trial -0.722 -0.038 0.031
Objct_cn:C_ -0.057 -0.071 0.045 0.041
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
The random effect structure seems to be too complex for the amount of
data. However, the variance for the (Intercept) under
NestID is 0 and the variance for
group_dummy under GroupID is 0.
This suggests that both nest as differences between groups contribute
minimally to the variance in the outcome model. Let’s drop both
effects.
zoi_model2 <- lmer(Zoi_duration ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
summary(enter_model2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data_enter
REML criterion at convergence: 4110
Scaled residuals:
Min 1Q Median 3Q Max
-5.7472 -0.0958 -0.0198 0.0611 19.4973
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 1.196e+02 10.93614
group_dummy 3.683e-04 0.01919 1.00
Residual 1.722e+02 13.12060
Number of obs: 506, groups: Bird_ID, 67
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.2873 1.2912 247.2062 4.095 5.72e-05 ***
Object_contrast -0.8899 1.1804 435.9864 -0.754 0.4513
Context_contrast -1.7996 1.7838 70.9217 -1.009 0.3165
Trial -0.5942 0.2595 438.1047 -2.289 0.0225 *
Object_contrast:Context_contrast 1.9818 2.3611 435.9877 0.839 0.4017
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_ Trial
Objct_cntrs 0.057
Cntxt_cntrs -0.433 -0.040
Trial -0.722 -0.038 0.031
Objct_cn:C_ -0.057 -0.071 0.045 0.041
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
check_model(enter_model2)
Seems like the interaction between object and context is non-significant, let’s drop it. Given the non-normal distribution, let’s logtransform the data as well:
summary(enter_model3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data
REML criterion at convergence: 696
Scaled residuals:
Min 1Q Median 3Q Max
-2.4872 -0.5604 -0.0780 0.4192 8.9279
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 0.083162 0.28838
group_dummy 0.009698 0.09848 1.00
Residual 0.196081 0.44281
Number of obs: 506, groups: Bird_ID, 67
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.143564 0.043897 223.101604 26.051 <2e-16 ***
Object_contrast 0.006807 0.039692 438.913282 0.171 0.864
Context_contrast -0.074761 0.046032 93.240157 -1.624 0.108
Trial -0.114456 0.008698 438.594691 -13.159 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_
Objct_cntrs 0.052
Cntxt_cntrs -0.324 -0.045
Trial -0.711 -0.035 0.038
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
There are still issues with the data distribution, let’s boxcox_transform:
Comparing the different models, the logtransfromed comes out as best:
anova(zoi_model, zoi_model2, zoi_model3, zoi_model4_boxcox)
refitting model(s) with ML (instead of REML)
Data: data
Models:
zoi_model3: log(Zoi_duration) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
zoi_model2: Zoi_duration ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
zoi_model4_boxcox: Latency_to_zoi_trans ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
zoi_model: Zoi_duration ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
zoi_model3 8 1342.4 1376.2 -663.21 1326.4
zoi_model2 9 6314.2 6352.3 -3148.11 6296.2 0.0 1 1
zoi_model4_boxcox 9 1416.4 1454.4 -699.18 1398.4 4897.9 0
zoi_model 11 6313.3 6359.8 -3145.63 6291.3 0.0 2 1
So we end up with these models:
# Time near object
summary(zoi_model3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Zoi_duration) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
Data: data
REML criterion at convergence: 1343.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.95388 -0.63743 -0.01154 0.62983 2.78860
Random effects:
Groups Name Variance Std.Dev. Corr
Bird_ID ind_dummy 0.2019 0.4493
group_dummy 0.2805 0.5296 1.00
Residual 0.6831 0.8265
Number of obs: 506, groups: Bird_ID, 67
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.38246 0.09114 174.24034 48.086 <2e-16 ***
Object_contrast 0.06530 0.07394 438.97006 0.883 0.378
Context_contrast 0.01085 0.07456 357.90017 0.145 0.884
Trial 0.02323 0.01618 438.49089 1.436 0.152
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Objct_ Cntxt_
Objct_cntrs 0.045
Cntxt_cntrs 0.036 -0.047
Trial -0.636 -0.035 0.042
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
This model considers the combined effect of latency and different contrasts:
latency_model <- lmer(log(Latency) ~ eat_vs_leave_contrast * Object_contrast * Context_contrast + Trial +
(- 1 + group_dummy | GroupID) +
(- 1 + ind_dummy + group_dummy + eat_vs_leave_contrast | Bird_ID),
data = data)
Compute and compare estimated marginal means for the
eat_model:
emm_eat <- emmeans(eat_model, ~ Object_contrast * Context_contrast)
pairs(emm_eat, adjust = "bonferroni")
Finally, a model is built for the Zoi_duration
variable:
zoi_model <- lmer(Zoi_duration ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(1 + group_dummy | GroupID) +
(1 + ind_dummy + group_dummy | Bird_ID),
data = data)
summary(zoi_model)
check_model(zoi_model)
As mentioned in the research report (RR), we may also consider analyzing latencies combined in a multivariate model:
multivariate_model
(Note: Details of the multivariate model should be added based on the specific analysis required.)
This document provides a comprehensive workflow for analyzing neophobia data, from data loading and preprocessing to statistical modeling and result interpretation. Further steps can be added as needed to refine and extend this analysis.